با استفاده از داده های زلزله ها در ایران و جهان به سوالات زیر پاسخ دهید.

توابع، داده ها و کتابخانه های زیر در طول تمرین مورد استفاده قرار می گیرند.

library(readr)
library(dplyr)
library(stringr)
library(ggmap)
library(highcharter)
library(plotly)
library(ggplot2)
library(rworldmap)
library(countrycode)
library(cowplot)
library(sp)
library(gganimate)
library(lubridate)

earthquakes = read_rds("data/historical_web_data_26112015.rds")

# convert coordinate to country
coords2country = function(points)
{  
  countriesSP = getMap(resolution='high')
  pointsSP = SpatialPoints(points, proj4string=CRS(proj4string(countriesSP)))  
  indices = over(pointsSP, countriesSP)
  indices$ISO3
}

۱. با استفاده از داده های historical_web_data_26112015.rds و استفاده از نمودار پراکنش سه بعدی بسته plotly نمودار طول، عرض و عمق زلزله ها را رسم نمایید. علاوه بر آن بزرگی هر نقطه را برابر بزرگی زمین لرزه قرار دهید.

plot_ly(earthquakes,
        x = ~Longitude, y = ~Latitude, z = ~Depth, size = ~Magnitude,
        marker = list(color = ~Magnitude,
                      symbol = 'circle',
                      sizemode = 'diameter',
                      colorscale = c('#FFE8A8', '#884848'),
                      showscale = TRUE),
        sizes = c(0.8, 8.8),
        text = ~paste('Province: ', Province,
                      '<br>City: ', City,
                      '<br>Earthquake Magnitude: ', Magnitude)) %>%
  add_markers() %>%
  layout(scene = list(title = 'Iran Earthquakes',
                      xaxis = list(title = 'Longitude', range = c(40, 80)),
                      yaxis = list(title = 'Latitude', range = c(20, 40)),
                      zaxis = list(title = 'Depth', range = c(0, 180))))

***

۲. پویانمایی سونامی های تاریخی را بر حسب شدت بر روی نقشه زمین رسم نمایید.(از داده زلزله های بزرگ استفاده نمایید.)

disaster = read_delim("data/disaster.txt",
                      "\t",
                      trim_ws = TRUE) %>%
  select(lon = LONGITUDE,
         lat = LATITUDE, 
         mag = INTENSITY,
         dep = FOCAL_DEPTH,
         country = COUNTRY,
         year = YEAR,
         t_death = TOTAL_DEATHS) %>% 
  na.omit()
eq_an = ggplot() +
  borders("world", colour="black", fill="gray") +
  geom_point(aes(x = disaster$lon,
                 y = disaster$lat,
                 size = disaster$mag,
                 frame = disaster$mag),
             color="blue") +
  ggtitle("Earthquakes In The World") +
  xlab("Longitute") +
  ylab("Latitude")

gganimate(eq_an, "eq_an.gif")

***

۳. نمودار چگالی دو بعدی زلزله های تاریخی ایران را رسم کنید.( از داده iran_earthquake.rds و لایه stat_density_2d استفاده نمایید).

eq_ir = read_rds("data/iran_earthquake.rds") %>% 
  filter(Long < 70)

ggplot(eq_ir, aes(x=Long, y=Lat)) +
  geom_point() +
  stat_density_2d(aes(fill = ..level.., colour = ..level..), geom = "polygon") +
  ggtitle("IR Earthquakes Density") +
  xlab("Longitude") +
  ylab("Latitude")


۴. احتمال اینکه در ایران در پنج سال آینده زلزله به بزرگی هفت ریشتر رخ دهد را محاسبه کنید. (از احتمال شرطی استفاده کنید.)

برای این سوال احتمال اینکه در پنج سال آینده زلزله رخ دهد به شرط آن که بیشتر از ۷ ریشتر باشد را محاسبه می کنیم. برای اینکار احتمال اینکه در پنج سال آینده یک زلزله رخ دهد را به شرط آنکه آن زلزله بزرگ تر از ۷ ریشتر باشد را بایستی محاسبه کنیم. در پنج سال آینده باتوجه به داده ها با احتمال ۱۰۰٪ زلزله رخ خواهد داد پس بایستی فقط احتمال اینکه در پنج سال آینده زلزله ی به بزرگی بیشتر از ۷ ریشتر رخ دهد را محاسبه کنیم.

eq_ir_7 = disaster %>% 
  filter(mag >= 7, country == 'IRAN') %>% 
  arrange(-year) %>% 
  mutate(dis = lag(year),
         dis = dis - year,
         dis = ifelse(is.na(dis),2018-year,dis))

eq_ir_7_under5 = eq_ir_7 %>% 
  mutate(uder5 = ifelse(dis <= 5, 1, 0)) %>% 
  group_by(uder5) %>% 
  summarise(count = n())

# +3 for continues < 5 years
prob = eq_ir_7_under5 %>% 
  summarise(prob =((7 + 3)/sum(count)))

knitr::kable(prob, 'html')
prob
0.625

۵. بر اساس داده های زلزله های بزرگ ابتدا تعداد و متوسط کشته زلزله ها را بر حسب کشور استخراج نمایید. سپس نمودار گرمایی تعداد کشته ها را بر روی کره زمین رسم نمایید.(مانند مثال زیر!)

disaster_country = disaster %>%
  group_by(country) %>%
  summarise(death_count = sum(as.integer(t_death)),death_mean = mean(as.integer(t_death))) %>%
  arrange(-death_count) %>% 
  mutate(country_code = countrycode(country,
                               "country.name",
                               "iso3c"))

knitr::kable(disaster_country %>% select(-country_code))
country death_count death_mean
CHINA 609347 18465.06061
IRAN 298903 15731.73684
ARMENIA 187000 62333.33333
ITALY 149711 7129.09524
PAKISTAN 146342 24390.33333
TURKMENISTAN 110001 55000.50000
PERU 83055 2373.00000
AZERBAIJAN 80087 26695.66667
ECUADOR 71202 17800.50000
CHILE 65469 2618.76000
TURKEY 59062 4218.71429
INDIA 55233 6904.12500
VENEZUELA 27808 3972.57143
GUATEMALA 23054 4610.80000
TAJIKISTAN 15529 5176.33333
MOROCCO 13728 6864.00000
RUSSIA 12025 1717.85714
NEPAL 11771 3923.66667
PHILIPPINES 10883 473.17391
MEXICO 10349 862.41667
ALGERIA 7308 1461.60000
JAPAN 6125 322.36842
GUADELOUPE 5000 5000.00000
UZBEKISTAN 4880 4880.00000
GREECE 3133 120.50000
YEMEN 2800 2800.00000
TAIWAN 2325 332.14286
INDONESIA 2202 100.09091
COLOMBIA 2158 269.75000
EL SALVADOR 1328 265.60000
MACEDONIA 1083 361.00000
AFGHANISTAN 934 155.66667
USA 655 27.29167
ARGENTINA 516 103.20000
KAZAKHSTAN 451 225.50000
GUINEA 443 443.00000
GEORGIA 280 93.33333
NEW ZEALAND 258 86.00000
ALBANIA 243 121.50000
PAPUA NEW GUINEA 180 30.00000
USA TERRITORY 168 84.00000
MONTENEGRO 131 131.00000
BULGARIA 130 65.00000
PORTUGAL 113 113.00000
COSTA RICA 93 31.00000
KYRGYZSTAN 79 39.50000
ETHIOPIA 70 35.00000
AZORES (PORTUGAL) 69 69.00000
SOLOMON ISLANDS 62 20.66667
MONGOLIA 30 30.00000
CANADA 28 28.00000
GHANA 22 22.00000
SLOVENIA 22 11.00000
IRAQ 20 20.00000
ROMANIA 18 6.00000
AUSTRALIA 12 12.00000
EGYPT 12 12.00000
UKRAINE 11 11.00000
DOMINICAN REPUBLIC 8 4.00000
CROATIA 5 2.50000
HONDURAS 3 3.00000
MYANMAR (BURMA) 3 3.00000
SERBIA 3 1.50000
BELGIUM 2 2.00000
CYPRUS 2 2.00000
DJIBOUTI 2 2.00000
AUSTRIA 1 1.00000
BOSNIA-HERZEGOVINA 1 1.00000
NETHERLANDS 1 1.00000
SOUTH AFRICA 1 1.00000
THAILAND 1 1.00000
mapCountryData(joinCountryData2Map(disaster_country, nameJoinColumn = "country_code"),
               nameColumnToPlot = "death_count", 
               mapTitle="Earthquakes Total Death")
70 codes from your data successfully matched countries in the map
1 codes from your data failed to match with a country code in the map
175 codes from the map weren't represented in your data

mapCountryData(joinCountryData2Map(disaster_country, nameJoinColumn = "country_code"),
               nameColumnToPlot = "death_mean", 
               mapTitle="Earthquakes Death Rate")
70 codes from your data successfully matched countries in the map
1 codes from your data failed to match with a country code in the map
175 codes from the map weren't represented in your data


۶. با استفاده از داده لرزه های بزرگ و به وسیله طول، عرض، شدت، عمق مدلی برای پیش بینی تعداد کشته های زلزله بیابید.

مدل glm را باتوجه به نوع داده ها با خانواده ی گاما ایجاد می کنیم و آن را ارزیابی می کنیم و درمی یابیم که به عمق و شدت آن بستگی دارد و مکان آن اهمیتی در مدل ندارد پس مدل جدید را براین مبنا مدل سازی می کنیم و ارزیابی می کنیم که نتیجه ی قابل قبولی به دست می دهد.

eq_glm = glm(t_death ~ mag + dep + lat + lon, family = Gamma(link = "inverse"), data = disaster)
summary(eq_glm)

Call:
glm(formula = t_death ~ mag + dep + lat + lon, family = Gamma(link = "inverse"), 
    data = disaster)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-4.055  -3.273  -2.611  -1.587   7.681  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.158e-03  2.561e-04   4.524 7.90e-06 ***
mag         -1.070e-04  2.406e-05  -4.445 1.12e-05 ***
dep          7.511e-06  3.933e-06   1.910   0.0568 .  
lat         -9.875e-07  2.824e-06  -0.350   0.7267    
lon         -2.387e-07  4.959e-07  -0.481   0.6304    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Gamma family taken to be 14.9034)

    Null deviance: 3965.6  on 429  degrees of freedom
Residual deviance: 3398.8  on 425  degrees of freedom
AIC: 5834

Number of Fisher Scoring iterations: 10
eq_glm = glm(t_death ~ mag + dep, family = Gamma(link = "inverse"), data = disaster)
summary(eq_glm)

Call:
glm(formula = t_death ~ mag + dep, family = Gamma(link = "inverse"), 
    data = disaster)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-3.895  -3.270  -2.619  -1.607   7.898  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.223e-03  2.326e-04   5.259 2.30e-07 ***
mag         -1.153e-04  2.197e-05  -5.247 2.44e-07 ***
dep          6.468e-06  1.967e-06   3.287   0.0011 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Gamma family taken to be 15.12215)

    Null deviance: 3965.6  on 429  degrees of freedom
Residual deviance: 3409.4  on 427  degrees of freedom
AIC: 5832.1

Number of Fisher Scoring iterations: 10

۷. با استفاده از داده worldwide.csv به چند سوال زیر پاسخ دهید. تحقیق کنید آیا می توان از پیش لرزه، زلزله اصلی را پیش بینی کرد؟

از روی داده ها پیش لرزه ها را استخراج کرده و روی آن ها تعداد و میانگین آن ها و قدرت زلزله ی اصلی را درکنار هم جمع کرده و مدلی را استخراج می کنیم که این پیش بینی را صورت دهد که خروجی آن نتیجه ی قابل قبولی را گزارش می کند.

worldwide = read_csv("data/worldwide.csv")

worldwide = worldwide %>%
  mutate(year = year(time),
         month = month(time), 
         day = day(time),
         location = as.data.frame(str_split_fixed(place, ',',2))$V2,
         location = ifelse(location == '' , as.character(place), as.character(location))) %>% 
  select(time, year, month, day, latitude, longitude, depth, mag, location) %>% 
  group_by(year, month, location) %>%
  arrange(time) %>% 
  mutate(mag_max = max(mag))
worldwide$groupid = worldwide %>% group_indices(year, month, location)

worldwide_temp = worldwide %>% 
  filter(mag == mag_max) %>% 
  ungroup() %>% 
  select(groupid,time_max = time)

worldwide = merge(worldwide, worldwide_temp)

worldwide = worldwide %>% 
  mutate(pre_past = ifelse(time > time_max, 'pre',ifelse(time == time_max, 'eq', 'past')))

worldwide_pred = worldwide %>% 
  filter(pre_past == 'pre' | pre_past == 'eq') %>% 
  ungroup() %>% 
  select(groupid, mag, pre_past) %>% 
  group_by(groupid,pre_past) %>% 
  summarise(avg_mag = mean(mag), count = n()) %>% 
  ungroup()

wwpm = worldwide_pred %>% 
  filter(pre_past == 'pre') %>% 
  select(groupid, count, avg_mag)

wwpmm = worldwide_pred %>% 
  filter(pre_past == 'eq') %>% 
  select(groupid, main = avg_mag)

x = merge(wwpm,wwpmm) %>% 
  select(-groupid)

glm_eq_pred = glm(count~., data = x)
summary(glm_eq_pred)

Call:
glm(formula = count ~ ., data = x)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
 -92.62   -15.38    -1.20    10.15  1100.56  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   14.463      6.533   2.214   0.0269 *  
avg_mag      -44.653      2.148 -20.790   <2e-16 ***
main          37.713      1.681  22.429   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 2502.604)

    Null deviance: 6924258  on 2218  degrees of freedom
Residual deviance: 5545771  on 2216  degrees of freedom
AIC: 23666

Number of Fisher Scoring iterations: 2

۸. گزاره " آیا شدت زلزله به عمق آن بستگی دارد" را تحقیق کنید؟ (طبیعتا از آزمون فرض باید استفاده کنید.)

cor.test(worldwide$depth, worldwide$mag, method = 'spearman')

    Spearman's rank correlation rho

data:  worldwide$depth and worldwide$mag
S = 4.176e+13, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.2376803 
mat_eq = as.matrix(worldwide %>%
                     select(depth, mag) %>% 
                     mutate(depth = depth + abs(min(depth)))
                   ) 

chisq.test(mat_eq)

    Pearson's Chi-squared test

data:  mat_eq
X-squared = 489570, df = 69011, p-value < 2.2e-16

تست وابستگی بیان می کند که غیرمرتبط بودن این دو متغیر از یکدیگر رد می شود پس می توان نتیجه گرفت که بین آن ها ارتباط وجود دارد و تست استقلال نیز بیان می کند که این دو متغیر نیز مستقل بودنشان از هم رد می شود و وابستگی آن ها نیز تایید می شود. پس فرض مستقل بودن این دو متغیر از هم رد می شود.


۹. میانگین سالانه زلزله ها را بر حسب کشور به دست آورید. آیا میتوان دلیلی در تایید یا رد تئوری هارپ ارائه کرد.

long_lat = worldwide %>%
  select(longitude, latitude)
long_lat$country = coords2country(long_lat)

worldwide$country_code = long_lat$country

country_code = worldwide %>% 
  ungroup() %>% 
  group_by(location) %>% 
  select(location, country_code) %>% 
  na.omit() %>% 
  distinct(country_code, location) %>% 
  ungroup() %>% 
  select(location, c_code = country_code) %>% 
  group_by(location) %>% 
  slice(1) %>% 
  ungroup()

worldwide_fixed_country_code = merge(worldwide, country_code)
worldwide_fixed_country_code = worldwide_fixed_country_code %>% 
  mutate(c_code = ifelse(!is.na(country_code),
                         as.character(country_code),
                         as.character(c_code)
                         ),
         country_code = c_code
         ) %>% 
  select(-c_code)

worldwide = worldwide_fixed_country_code

worldwide_accu = worldwide %>%
  group_by(year, country_code) %>% 
  summarise(eq_count = n(),
            mean_mag = mean(mag),
            mean_dep = mean(depth))

eq_2015 = worldwide_accu %>%
  filter(year == 2015) %>% 
  select(eq_count)

eq_2016 = worldwide_accu %>%
  filter(year == 2016) %>% 
  select(eq_count)
                    
eq_2017 = worldwide_accu %>%
  filter(year == 2017) %>% 
  select(eq_count)
                   
eq_2018 = worldwide_accu %>%
  filter(year == 2018) %>% 
  select(eq_count)

t.test(eq_2015$eq_count,eq_2016$eq_count)

    Welch Two Sample t-test

data:  eq_2015$eq_count and eq_2016$eq_count
t = -2.462, df = 109.01, p-value = 0.01538
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -381.34944  -41.19174
sample estimates:
mean of x mean of y 
 32.98413 244.25472 
t.test(eq_2015$eq_count,eq_2017$eq_count)

    Welch Two Sample t-test

data:  eq_2015$eq_count and eq_2017$eq_count
t = -2.1236, df = 108.3, p-value = 0.03598
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -387.70921  -13.36027
sample estimates:
mean of x mean of y 
 32.98413 233.51887 
t.test(eq_2015$eq_count,eq_2018$eq_count)

    Welch Two Sample t-test

data:  eq_2015$eq_count and eq_2018$eq_count
t = -1.2653, df = 96.949, p-value = 0.2088
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -214.29194   47.43411
sample estimates:
mean of x mean of y 
 32.98413 116.41304 

با توجه به اینکه پروژه ی هارپ در سال ۲۰۱۵ متوقف شد زلزله های این سال را با سال های بعدی آن مقایسه کردم و تست t را بر روی آن ها اجرا کردم تا تفاوت این سال ها را باهم مقایسه کنم. نتایج این تست نشان داد که این سال ها بایکدیگر تفاوت معناداری ندارند و این تئوری را می توان رد کرد.


۱۰. سه حقیقت جالب در مورد زلزله بیابید.

۱۰.۱. نمودار جعبه ای زلزله ها در سال ها و جایگاه ایران در آن

ww_mag_year = worldwide %>% 
  group_by(year, country_code) %>% 
  summarise(mean_mag = mean(mag)) %>% 
  select(year, mean_mag) %>% 
  ungroup()

iran_w = worldwide %>% 
  group_by(year, country_code) %>% 
  summarise(mean_mag = mean(mag)) %>% 
  filter(country_code == 'IRN') %>% 
  select(year, mean_mag) %>% 
  ungroup()
  
ggplot(ww_mag_year,aes(x = year, y = mean_mag, group = year)) +
  geom_boxplot() +
  geom_line(data = iran_w,aes(x = year, y = mean_mag, group = 1, color = 'Iran'), size = 1.5) + 
  ggtitle("Average Earthquake Magnitude In Worldwide") + 
  ylab("Magnitude Mean") +
  xlab("Year")

۱۰.۲. ده زلزله خیزترین کشورها

country_reg = country_code %>% 
  group_by(c_code) %>% 
  summarise(reg_num = n()) %>% 
  ungroup() %>% 
  select(country_code = c_code, reg_num)

ww_mag_count = worldwide %>% 
  group_by(country_code) %>% 
  summarise(mean_mag = mean(mag), count = n()) 

ww_mag_count = merge(ww_mag_count, country_reg)

top_10 = ww_mag_count %>% 
  mutate(score = (count*mean_mag)/reg_num) %>% 
  arrange(-score) %>% 
  slice(1:10)

knitr::kable(top_10, "html")
country_code mean_mag count reg_num score
JPN 4.531295 2486 1 11264.80
PNG 4.576130 2191 1 10026.30
NZL 4.537366 2058 1 9337.90
PRI 2.886450 2972 1 8578.53
PHL 4.531641 1438 1 6516.50
IDN 4.462053 4277 3 6361.40
TON 4.628141 1329 1 6150.80
VGB 2.986832 1632 1 4874.51
VUT 4.580264 1059 1 4850.50
RUS 4.425916 1092 1 4833.10

۱۰.۳. نحسی عدد ۱۳ در وقوع زلزله در روز ۱۳ام.

thirteen = worldwide %>% 
  filter(day == 13)

not_thirteen = worldwide %>% 
  filter(day != 13)

t.test(thirteen$mag, not_thirteen$mag, alternative = "greater")

    Welch Two Sample t-test

data:  thirteen$mag and not_thirteen$mag
t = 6.5195, df = 2091, p-value = 4.407e-11
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
 0.0961231       Inf
sample estimates:
mean of x mean of y 
 3.873674  3.745097 

این تست رد می شود و نمی توان نتیجه گرفت که روز ۱۳ام روز نحسی برای زلزله است.